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Abstract —In this paper, a sparse-based method for the es¬ 
timation of the parameters of multidimensional (i?-D) modal 
(harmonic or damped) complex signals in noise is presented. The 
problem is formulated as R simultaneous sparse approximations 
of multiple 1-D signals. To get a method able to handle large 
size signals while maintaining a sufficient resolution, a multigrid 
dictionary refinement technique is associated with the simulta¬ 
neous sparse approximation problem. The refinement procedure 
is proved to converge in the single R-D mode case. Then, for 
the general multiple modes R-D case, the signal tensor model 
is decomposed in order to handle each mode separately in an 
iterative scheme. The proposed method does not require an 
association step since the estimated modes are automatically 
“paired”. We also derive the Cramer-Rao lower bounds of the 
parameters of modal R-D signals. The expressions are given in 
compact form in the single R-Y> mode case. Finally, numerical 
simulations are conducted to demonstrate the effectiveness of the 
proposed method. 

Index Terms —Multidimensional modal retrieval, frequency 
estimation, simultaneous sparse approximation, multigrid dictio¬ 
nary refinement, Cramer-Rao lower bound, harmonic retrieval 


I. Introduction 

T he problem of estimating the parameters of sinusoidal 
signals from noisy measurements is an important topic in 
signal processing and several parametric and nonparameteric 
approaches have been developed for one-dimensional (1-D) 
signals [1]. Recently, this problem has received a renewed 
interest thanks to the emergence of multidimensional (i?-D) 
applications. Indeed, parameter estimation from R-D signals 
is required in numerous applications in signal processing and 
communications such as nuclear magnetic resonance (NMR) 
spectroscopy, wireless communication channel estimation [2] 
and MIMO radar imaging [3]. In all these applications, signals 
are assumed to be a superposition of R-D sinusoids or, more 
generally, of exponentially decaying i?-D complex exponen¬ 
tials (modal signals). As for the 1-D case, the crucial step is 
the estimation of the R-D modes (including frequencies and 
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damping factors) because they are nonlinear functions of the 
data. 

In order to achieve high resolution estimates, parametric 
approaches are often preferred to nonparamatric ones. Several 
parametric R-D methods (i? > 2) have been proposed. 
They include linear prediction-based methods such as 2- 
D TLS-Prony [4], and subspace approaches such as matrix 
enhancement and matrix pencil (MEMP) [5], 2-D ESPRIT [6], 
multidimensional folding (MDE) [7], improved multidimen¬ 
sional folding (IMDE) [8], [9], Tensor-ESPRIT [10], principal- 
singular-vector utilization for modal analysis (PUMA) [11], 
[12] and the methods proposed in [13], [14]. All these methods 
perform at various degrees but it is generally admitted that 
they yield accurate estimates at high SNR scenarios and/or 
when the frequencies are well separated. This is obtained at 
the expense of computational effort. Eor instance, MDE, IMDE, 
2-D ESPRIT and MEMP are ESPRIT-type techniques. They 
require to build large size matrices and apply the ESPRIT- 
based method, which make their computational complexity 
very high particularly in the case of large R-D signals. The 
Tensor-ESPRIT algorithm uses the structure inherent in the 
i?-D data at the expense of a high computational complexity. 
Recently, TPUMA [11] was proposed as an accurate and 
computationally efficient multidimensional harmonic retrieval 
method, which attains the Cramer-Rao lower bound (CRLB) 
and does not require to build large size matrix or tensor. 
However its performance degrades rapidly with the increase 
of the number of components present in the R-D signal. 

Recently, methods based on sparse approximations have 
been proposed to address the harmonic or modal retrieval prob¬ 
lem [15]-[22]. Eor time-data spectral estimation, the dictionary 
is formed from a set of (normalized) complex exponentials 
potentially embedded in the data, which allows one to easily 
include some prior knowledge about the position of certain 
known modes. More generally, the usual choice is a uniform 
spectral grid obtained by sampling the frequency and damping 
factor lines. Clearly, a fine grid will lead to a good resolution 
but, on the other hand, it will result in a huge dictionary [15]. 
This complexity is further increased in the case of R-D signals 
in which we are confronted with 2R-D grids. 

In order to reduce the computational burden, a multigrid 
scheme for sparse approximation was proposed in [19] to 
iteratively refine the dictionary starting from a coarse one. At 
each iteration, a sparse approximation is performed and then 
new grid points (atoms) are inserted in the vicinity of active 
ones leading to a multiresolution-like scheme. So the multigrid 
algorithm in [19] refines jointly R 2-D grids. This algorithm 
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is efficient but has mainly two drawbacks: 1) it does not have 
convergence guarantees, 2) the dictionary becomes intractable 
for large signals when R> 2. 

The goal of the present paper is to propose a fast multi¬ 
dimensional modal estimation technique able to handle large 
signals and yielding a good estimation accuracy. 

1) First, the proposed approach, as for some parametric 
methods for modal retrieval, is based on the idea of esti¬ 
mating the parameters independently along each dimension 
r = 1,..., i?. It will be shown that the simultaneous sparse 
approximation concept [20], [23] is well-suited for R-D 
modal retrieval (R > 2). 

2) The second contribution consists in the proposition of a 
new multigrid scheme which amounts to consider a two- 
step refinement of 1-D grids, the first step for frequencies 
and the second one for damping factors. One advantage of 
the two-step multigrid is that it reduces the computational 
time. The convergence issue of the proposed multigrid 
strategy is analyzed firstly for single tone {F = 1) case, 
and convergence conditions are derived. Condition for 
convergence are expressed in terms of atom positions in 
the initial dictionaries. 

3) The extension of this result to the multiple tones case 
(F > 1) is not trivial because, not only it depends on 
the selected sparse approximation algorithm, but also on 
the coherence of the dictionary [23]. Indeed, due to the 
refinement strategy, the resulting dictionary is far from 
being uncorrelated which may prevent convergence even 
in the noiseless case. Consequently, for the multiple tones 
case, we exploit an alternative representation of the data 
model enabling the extraction of the R-D signal tones 
separately. Therefore, the third contribution of this paper 
consists in deriving a new algorithm for estimating param¬ 
eters of i?-D damped signals in which the results of the 
previous contribution apply. The effectiveness of the new 
algorithm for multiple R-D tones is also analyzed. One 
very interesting by-product of this approach is that the 
pairing of R-D parameters is achieved for free, without 
any further association stage. 

To assess the performances of an estimation method, the 
usual way consists in comparing the variance of the estimates 
to the CRLB. In [5] Y. Hua derived the CRLB for 2-D 
frequencies, i.e., undamped 2-D exponentials; no damped 
signals are considered. Closed-form expressions of the CRLB 
for the general undamped R-D case are derived in [24]. CRLB 
for 2-D damped signals are derived in [25]. Therefore, to the 
best of our knowledge, no compact expressions of the CRLB’s 
are available for the general R-D modal (damped) signal. 
Thus, another contribution of this paper is the derivation of 
the CRLB’s for the frequency, damping factor, amplitude and 
phase of the R-D modal signal. 

The remainder of this paper is organized as follows. In 
section II, we introduce notation and present the R-D modal 
retrieval problem. In section III, we formulate the R-D modal 
estimation problem as R simultaneous sparse estimation prob¬ 
lems, show how to construct a modal dictionary on a uniform 
grid and then describe the new fast multigrid strategy. In 


section IV, we give sufficient conditions for convergence of the 
multigrid dictionary refinement scheme in the case of single 
tone i?-D signals. In light of these new results, we propose 
in section V a new efficient algorithm for multiple tones R- 
D modal signals. In section VI, we derive the expressions of 
the CRLB’s for the parameters of R-D damped exponentials 
in Gaussian white noise. We then give the CRLB in the 
cases of single damped and undamped R-D cisoids. The 
effectiveness of the proposed method is demonstrated using 
simulation signals in section VII. Finally, conclusions are 
drawn in section VIII. 

IT Notation and Problem Statement 

A. Notation 

In this paper, scalars are denoted as lower-case letters 
(a, 5, a), column vectors as lower-case bold-face letters (a, b), 
matrices as bold-face capitals (A, B), and tensors as calli¬ 
graphic bold-face letters (^, B). Let (•)^, (•)'^ and (•)^ denote 
the transpose, the Hermitian transpose and the pseudo-inverse, 
respectively. The symbols “0” and “Kl” will denote the Khatri- 
Rao product (column-wise Kronecker) and the Kronecker 
product, respectively. Both words “mode” and “tone” are used 
to refer to a component of the multidimensional signal. The 
tensor operations used here are consistent with [26]: 

• the outer product of two tensors A. G and 

B G is given by: 

C = A0 B G C^^'^-xMRXKiX-xKff ^ 
c(mi,...,m_R,A:i,...,fcjv) = 

a{mi,...,mR)b{ki,...,kN) (1) 

• the contraction product acting on the rth index of a tensor 

A G and the 2nd index of a matrix U G 

^KxMr is; 

r ^ 

. . . ,mr-l,kr,mr+l, . . -,771^) = 

Mr 

■. ,mR)u{kr,mr) (2) 

nir — 'i 

• the matrix A(r) G £,Mrx(Mi---Mr-iMr+i---MR) j-gpresents 
the unfolding (dimension-r matricization) of the tensor A 
and corresponds to the arrangement of the dimension-r 
fibers of A to be the columns of the resulting matrix. 

• II A.|p denotes the the Frobenius norm for tensors: = 

|a(mi,...,mfl)|2. 

B. Problem Formulation 

An R-D modal signal is modeled as the superposition of F 
multidimensional damped complex sinusoids: 

F R 

y{mi, tur) = X! n + e{7ni,rriR) (3) 

/=! r=l 

where rrir = 1,..., Mr for r = 1,..., i?. Mr denotes the sam¬ 
ple support of the rth dimension, aj^r = exp (a/,r + j^f,r) G 
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C is the /th mode in the rth dimension, r=i’ 

a/,r G are the damping factors, {w/,r = ^=1 

are the angular frequencies, and c/ = Xfexp{j(j)f) is the 
complex amplitude of the /th mode where A/ = |c/| de¬ 
notes the magnitude and // the phase. e{mi,m 2 , ■ ■ ■ ,mu) 
is a zero-mean complex Gaussian white noise with variance 
cr^ and mutually independent components in all dimensions. 
Throughout this paper, the tilde symbol (') denotes a noisy 
signal; e.g. y{-) = y{-) -f e(-). 

In a tensor form, the R-D signal in (3) may be written as 

y = y + e ( 4 ) 

where {y,y,£} G £;MixM 2 x .xMr_ problem consists 
in estimating the set of parameters {a f,r}f’J^i ^=1 {cf}J^-^ 

from the R-D signal samples. 

III. Simultaneous Sparse Approximation eor R-D 
Modal Signals 

A. Tensor Formulation of the Data Model 

The noise-free data tensor y in (4) can be written in the 
following form: 

F 

y = ^ c/ a /4 (g) a /,2 G) • • • G) a/,_R (5) 

/=! 

where a/,^. = [1,a/,r, • ■.,, r = 1,... ,R. Equation 
(5) is called the Canonical Polyadic (CP) decomposition form, 
or the Candecomp/Parafac decompostion of the tensor y [26], 
[27]. The CP model (5) can be concisely denoted by 

= [c;Ai,A2,...,Ajil (6) 

where Ar = [ai^r, a 2 ,r, ■ ■ ■, au.r], r = 1, ■ ■ ■ ,R, and c = 
[ci, C 2 ,..., cfY is ths vector of complex amplitudes. Using 
these definitions, the matricized form of y along the rth 
dimension is given by 

= Ar Ac(A/j 0 • • • 0 Ar+l 0 Ar_l 0 • • • 0 Ai)"'' 

(7) 

where Ac = diag(c). Then, we can write 

Y(r) = ArHr -f E(r) (8) 

where G is 


considered as multiple experiences involving the same one¬ 
dimensional signal generated by the modes a/,r) / = 1) • ■ •) E, 
but with different amplitudes for each experience. This prop¬ 
erty will be used in the next section to formulate the problem 
of estimating the mode coordinates in the rth dimension as a 
simultaneous sparse approximation problem. 

B. Simultaneous Sparse Approximation 

Assuming' that Mr > F,\/r, it is easy to see from (10) 
that for a fixed r the mode coordinates {a/,r}/Li < F) 

in the rth dimension are identifiable from any column of 
Y(r)- This process can also be repeated on each dimension 
r = 1,..., i? to get all the modes coordinates. In practice, 
we have to replace the matrix Yj-^.) by its noisy counterpart 
Y(j.) accounting for the additive white noise. In this case, 
(10) holds only approximately. Consequently, for each column 
y(r),m',,; iii-r = 1) ■ • ■; 'he modal estimation problem can 
be formulated as a sparse approximation problem correspond¬ 
ing to the following constrained optimization: 

x„j; = argmin||x||o subject to - Qrx ||2 < e 

(11) 

where G N ^ Mr, is a (known) modal dictio¬ 

nary, X G is a (sparse) vector containing the coefficients of 
the activated columns in Q^., and e is a small reconstruction 
error related to the noise variance. The pseudo-norm ||x||o 
counts the number of nonzero elements in a vector x. The 
design of Q^. is discussed in section III-C. The fact that each 
vector y(r),m' corresponds to a 1-D signal generated by the 
same modes implies that the position of nonzero entries in 
Xm/ should be the same for mj, = 1,2,..., Mjj. Let X be 
the matrix defined by 

X = [xi,X 2 ,... ,xm;], (12) 

then the sparsity of X may be measured by computing the 
Euclidian norms of the rows; those providing a nonzero norm 
define the rows of the activated atoms (which are estimations 
of modes a/ ,, in the dimension r) in the dictionary Q^. 
Therefore, we are facing a simultaneous sparse approximation 
problem: 

Xr = argmin ||X|| 2 ,o subject to ||Y(,r) ~ QrX|||n < e 

(13) 

where 


= Ac(Afi 0 • • • 0 A^+i 0 Ar-i 0 • • • 0 Ai)"' (9) 

R 

and M^ = M^. Therefore 

k=l 

k^r 

"^(r) — [y{r),li ■ ■ • •iY{r),M^] 

F F 

= '^hr{f,l)af^r,- ■ ■,'^hr{f,Mr)af^r ( 10 ) 

_/=l /=! 

where hr{f,rn'f) is the {f,m'r) entry of the matrix H^, for 
/ = 1,..., U and m'r = 1,..., Mf We observe that, for a 
given r, the columns y(r),m' of Y(r) are linear combinations 
of the vectors {af^r}f-i- Hence, the columns y(r),m; can be 


||Y(,) - Q.XIII = ||vec(Y(,) - Q.X)||2, (14) 

||X||2.o= [||XTj|2 ••• ||XY:||2]T , ( 15 ) 

0 

and X"’- stands for the nth row of X. The simultaneous 
sparse representation models, called also Multiple Measure¬ 
ment Vectors (MMV), have been studied from several angles 
of view, and different approaches have been proposed [28], 
using greedy strategies [29] such as Simultaneous Orthogo¬ 
nal Matching Pursuit (SOMP) [23], convex relaxation meth¬ 
ods [30], randomized algorithms such as REduce MMV and 
BOost (ReMBo) [31] and subspace-augmented MUSIC [32]. 
As the goal of the present paper is to develop a fast approach 

'Note that this assumption is considered only in this section. 



4 


well adapted to large signals, we restrict our attention to the 
SOMP algorithm [23], reported in Appendix A. However, 
it is worth mentioning that, in more intricate cases and/or 
small size signals, much more efficient simultaneous sparse 
algorithms may be used at the price of an increased com¬ 
putational burden. A straightforward way to get the i?-tuples 
..., consists in estimating the modes a/,r in 

the R dimensions using matrices 'Y = 1,..., i?, which 
requires a further pairing step to form the i?-D modes in the 
multiple tones case {F > 1). To get accurate estimates using 
the described scheme, two conditions have to be satisfied, 1) 
the dictionary should contain all possible modes present in 
the signal, 2) the sparse approximation method should have 
sufficient guarantees for selecting the true atoms from the 
dictionary, which is known as “exact recovery guarantees”. 
These problems are discussed in the following sections and 
an alternative representation of the data is used to avoid the 
pairing stage in the multiple tones case. 

C. Modal Dictionary Design and Multigrid Strategy 

1) Uniform Modal Dictionary: The dictionary S 

^MrxN defined as follows. Let be the number of 

points of a uniform grid covering the frequency interval [0,1). 
Similarly, let Np be the number of points of a uniform grid 
covering the damping factor interval (/3min;0], where /^min is 
a lower bound on {ct/Then is given by 

Qr = [qr(0,0),... ,qr((JV^ - l)S^,0),qr(0,Sp),..., 

q,((N^ - 1)S^, 5p ),..., q,((iV^ - 1)S^, (Np - 1)5^)] 

( 16 ) 

where qr{fJ.,(3) = ar(/r,/3)/||ar(/r,/3)||2 with ar{ti, P) = 

[l^g(/3+j27rA.)^_^_^g(/3+i27rA.)(M.-l)]T, ^ P^i^/Np, and 

5^ = 1/A)i. In short, Q^. is obtained from a discretization 
of the (z4, a) plane. Each point of the grid corresponds to 
a hypothetic mode. The total number of columns in Qr is 
N = NfiNp ^ F, each of them is called atom. In the aim of 
reducing the computational complexity, we propose to estimate 
frequencies and then damping factors by calling twice the 
sparse approximation method. At the first step, the frequencies 
are estimated using a harmonic dictionary. In the second step, 
the damping factors are estimated using a modal dictionary 
formed by the already estimated frequencies and a damping 
factor grid. These two steps are explained in section IV. 

2) Multi-Grid Dictionary Refinement: To achieve a high- 
resolution modal estimation, a possible way is to define 
uniform grids as before and selecting very small values for 5^ 
and 5p to retrieve the frequencies and damping factors, respec¬ 
tively. As a consequence, the resulting dictionaries will lead to 
prohibitive calculation cost and memory capacities requested. 
Rather, we propose to start with a coarse one {N^ and Np 
low) and to adaptively refine it through a multigrid scheme. 
The procedure is the same for estimating the frequency and 
damping factor. The principle is sketched on Figure 1. The 
main idea is the adaptation of the dictionary as a function 
of the previous dictionary and the estimated coefficients. Let 
t be the current grid level (£ = 0,..., L — 1). At level £, 
we first restore the signal Xr(f) related to the dictionary 


Level 



Fig. 1. The multigrid dictionary refinement procedure with 77 = 1. (Q) atoms 
in the dictionary; ((•)) activated atoms; ((g)) new atoms 


Algorithm 1; Dictionary refinement (DICREF) 

input : A vector d G of sorted frequencies or damping factors, an 
index set H of activated atoms, the number of atoms 77 € N to 
add at each side of an activated one 
output: Updated vector dupdated 

for f = 1 : numel{Q) do 

dj,i = linspace (d(r2(2) — 1 ), 77 ) 

di 2 = linspace(d(r2(2)), d(r2(2) + 1 ), 77 ) 
di = : r])f 

end 

^updated “ Union(dl, . . . , ) 

return d^p^j^ted 


Qr{£) by applying the SOMP method. Then we refine the 
dictionary by inserting atoms inbetween pairs of Qr(f), in the 
neighborhood of each activated atom and we apply again the 
SOMP method at level £ -f 1 to restore Xr(£ -f 1) with respect 
to the refined dictionary Qr(£ +1). This process is repeated 
until a desired level of resolution is reached. Algorithm 1 
presents the one-step dictionary refinement (DICREF), from 
level £ to £ -f 1, where, for a and b reals, linspace(a, b, rf) 
generates a set of rj equispaced points in the interval [a, b]. The 
difference between the present framework and that in [19] is 
the following. In [19] the multigrid algorithm refines jointly 
R 2-D grids, which leads to expensive computations when 
R>2, without convergence guarantees. The present mutigrid 
scheme refines linear grids, which leads to low computational 
complexity with convergence guarantees as we will show in 
the next section. 

Finding the convergence conditions of the new multigrid 
strategy in the general case (multiple tones) is not easy and 
depends on the selected sparse approximation algorithm. By 
contrast, it is possible to show that, under mild conditions, the 
convergence may be guaranteed in the single tone case. This 
issue is discussed in the next section. In section V, we make 
use of an alternative representation of the data model in the 
case of multiple tones and a method allowing one to retrieve 
the signal tones separately. 

IV. Single R-D Mode Estimation 

In the previous section, we have shown how the R-D modal 
retrieval problem may be tackled using a sparse approximation 
algorithm by estimating the set of parameters in each dimen- 
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sion r = 1,. .., i?. Here, we give the sufficient conditions for 
convergence of the multigrid dictionary refinement scheme for 
F = 1. Without loss of generality, we set i? = 1. For notation 
simplicity, we omit reference to the dimension index r. 

According to (3), the 1-D modal signal containing a single 
mode can be written as follows; 

(17) 

Let Q be a normalized modal dictionary Q = [qi,..., qjv], 
with 

. 

qn = exp(/ 3 „ + G [ 0 , !),/?„ G (/ 3 min, 0 ], for 

n = The single tone sparse approximation of y 

with respect to Q is the solution of the criterion: 

min J(x) = ||y - Qx|p s.t. ||x||o = l. (19) 

X 

The optimal solution is given by 

< = ql^y, jv}\n = 0: = ||y||^-y'^qnq^y 

( 20 ) 

where n is the selected column number in Q. Finally, the 
minimum J(x*) is reached for an atom q„ that maximizes 
= y^qnqt^y = Iqt^yp, n = i,...,N. 


A. Estimating the Frequency: The Harmonic Dictionary 

First, we estimate frequency vi using a harmonic dictionary 
(i.e. assuming /3„ = 0, Vn). In this case, we have; 

J'iFn) = ■ '^ 21 ) 

The following theorem gives a sufficient condition for the 
multigrid dictionary refinement scheme to converge to the 
global maximum of J'. 

Theorem 1: Let y{m) be a single tone {F = 1) noiseless 
signal of length M and Q(f = 0) = [qi q 2 ... qAr(o)]^ be the 
initial harmonic dictionary in which the columns are sorted in 
increasing order of /i„(0),n = 1, 2,..., A^(0) and covering 
the frequency interval [0,1): /ii(0) = 0 and /rjv(o)(0) = 
1 — 1/M. Then the refinement scheme is convergent (i.e. 
3n G {1,..., N{i)} s.t. lim£_).oo = vi) if the following 
condition is satisfied: 

max |m„+i( 0) - Ain(0)| < 2 Cm (22) 
ne{l,....Af(0)-l} 

where C,m is a constant depending only on M. 

Proof It is easy to check that the global maximum of J'{iin) 
is reached for fin = ui, Vai. Figure 2 shows the variation of 
J'iUn) as a function of for i/i = 0.1, /3„ = 0 and M = 10. 
For Q!i = 0, J'{gLn) reduces to a Fejer kernel which has exactly 
one local maximum in the interval [z/i + fc/M, i/i + (fc + l)/M], 
k 0. Let J[ be the maximum value of J'{fj.n) in the interval 
[vi + 1/M, ui + 2/M] and i/i + (/m be the value of such 
that J'{iJLn = vi+(/m) = J'l in the interval [vi, vi + l/M] (we 



Fig. 2. in the single mode case with ui =0.1 and ^ 


assume^ that M > 2). For the dictionary refinement strategy to 
converge to the global maximum, it is sufficient to the sparse 
approximation algorithm to select, at a given level £, an atom 
whose frequency satisfies |p„*(f) — vi \ < (m < 1/lbf, where 
p„. {£) = argmax„ J'{nn)- Indeed, if p„. {£) G + 

(/m) then adding two atoms whose frequencies are located on 
both sides of /i„. {£) will lead to the selection, at level f + 1, 
of an atom that satisfies |/i„. (f +1) — z/i| < |p„. (.f) — Z 2 i|; the 
distance between the selected atom and the true frequency is a 
monotonically decreasing sequence. Finally, the convergence 
is guaranteed if the initial dictionary contains an atom n such 
that |pn(0) — i^i\ < Cm, which is satisfied if 

max |m™+i( 0) - Pn(0)| < 2Cm- (23) 
ne{i,....w(o)-i} 

given the fact that the sequence {/ira(O)} covers the interval 
[0,1). For ai < 0, the main lobe of J'{y,n) becomes broader 
and Cm larger than for ax = 0. Consequently, condition (23) 
is also sufficient for ai < 0. ■ 

Corollary 1: In the single tone case, the harmonic dictionary 
refinement is convergent if the initial frequency grid {£ = 0) 
is the Fourier grid. 

Proof Fourier bins are obtained for N = M and iin{0) = 
{n — 1)/M. Since Cm > 1/2M, the proof is straightforward 
because |Ain+i(0) - M„(0)| = 1/M < 2Cm- ■ 

It is important to note that condition (23) is sufficient but 
not necessary. Moreover, this condition is established when 
adding a single atom on both sides of the selected one (i.e. 
p = 1 in Algorithm 1). When p ^ 1, the condition may be 
relaxed and the rate of convergence is expected to be higher. 


B. Estimating the Damping Factor: The Modal Dictionary 


Assume that the previous sparse approximation method 
using a harmonic dictionary has converged to select an atom 
with /r„ = vi. Now, we have to estimate the damping factor 
ai- We form a modal dictionary using the damping factor 
grid and the frequency z^i, i.e. gn = exp(/3„ + j2'KVi). 
Consequently, 


J\Pn) 


|cip(l — /I ~ 

\ — g2/3TiA^ J 


(24) 


^The case of M < 2 in not of practical interest but the theorem is still valid 
by setting \ because J'{^n) is a monotonically decreasing function 

in the interval [t^i, minji, ^ + ^i}]- 
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Theorem 2: Let y{m) be a single tone {F = 1) noiseless 
signal of length M and Q(0) = [qi q 2 • • • qAf(o)]^ be 
the initial modal dictionary formed using the frequency 
i.e., Qn = exp(/3„(0) + j2Tri'i), where i^i is the frequency 
of signal y. The columns are sorted in increasing order of 
/3„(0),n = 1,2,...,TV and covering the damping factor 
interval (/3mitn0]. Then the refinement scheme is convergent 
(i.e. 3n s.t. limr^oo/3nW = cti) if cti G (/3min,0]. 

Proof Let g^Pn) be the derivative of J'{Pn) in (24) with 
respect to /3„. It is easy to check that g{Pn) > 0 for /3„ < cti, 
g{Pn) < 0 for /3„ > cti, and p(/3„) = 0 when /3„ = cti. 
In other words, J'{Pn) is monotonically increasing before 
the maximum reached at cti and monotonically decreasing 
after cti. Therefore, the multigrid algorithm converges to ai 
if Pmin ^ CXl- H 

As a consequence of Theorem 2, the initial modal dictionary 
can be formed using only two points in the damping factor 
grid: /3i(0) = /3i„in and ^ 2 ( 0 ) = 0. 

We can now state that the multigrid algorithm based on two 
sparse approximations (for frequency and then damping factor) 
converges in the single tone case under some conditions. Note 
that in the noisy case when the SNR is sufficiently high, the 
convergence analysis is still valid as in the noiseless case, 
and the proposed multigrid sparse scheme for single tone 
converges to the global maximum of the Fejer kernel. The 
extension to the single tone R-D modal retrieval problem 
is straightforward and can be performed according to the 
formulation presented in Section III-B. The details of this 
approach (STSM: Single Tone Sparse Method) are presented in 
Algorithm 2. The algorithm takes as input a noisy single tone 
R-D signal, and a couple of integers 77 ^ and r]a that correspond 
respectively to the number of frequency and damping factor 
atoms to be added on both sides of the corresponding selected 
ones. Next, for each dimension r = 1,..., i?, we execute two 
tasks to estimate the frequency and then the damping factor: 
in each step we apply SOMP combined to DICREF algorithm 
using corresponding dictionaries and taking into account the 
convergence conditions discussed previously. Then parameters 
of Or, i.e., I'r and ar, are given by the corresponding selected 
atoms. This approach will be exploited in the next section for 
the multiple tones case. 

V. Multiple R-D Modes Estimation 

In the multiple tones case, sparse approximation algorithms 
yield suboptimal solutions when the coherence of the dictio¬ 
nary is high [29]. This is a crucial point because the refinement 
procedure will increase the coherence with increasing £, which 
may prevent convergence even in the noiseless case. In the 
following, we present a low complexity algorithm that is ac¬ 
curate and robust in the presence of noise. The idea is to begin 
by an initialization step where F single tone modal signals 
of order i? — 1 are extracted from the R-D signal. Then an 
iterative technique is proposed to improve this decomposition 
and estimate the parameters. 

It is assumed that the frequencies are distinct in at least one 
dimension with Mr > F. Then dimensions are permuted such 


Algorithm 2: Sparse multigrid method for single tone 
estimation (STSM) 

input : A tensor y S ^ ^ , » 7 q, ) g N X N 
output: Parameters of the single R-D mode: ai,... ,a/^ 

initialization: {ki^,kct) = (0,0) 
initialize and using 

for r = 1 : i? do 

while halting criterion false do 

kjy = ki/ 1 

= SOMP{Q{dl''‘'\o), Y(,,),Iter = 1) 

^ DICREF(d],''''\q],'‘\»7,.) 

end 

while halting criterion false do 
ka — ka + 1 

= SOMP (q , Y(,)) 

d^^+i) = DICREF(dL''“\qi'‘\r;„) 

end 

end 

return ai ,..., an 


that the dimension with distinct frequencies becomes the first 
one (r = 1). 

A. From Multiple Tones to Multiple Single-Tone Signals 
According to (5), y can be written as 

y - Xr+i^F • Ar • (25) 

r=l i?+l 

= <S*Ai (26) 

where Xr^i^f is the diagonal tensor of order i? -f 1 and size 
F X F X ■ ■ ■ X F, containing ones on its diagonal, and 

S=Xr+i,f • Ar • (27) 

r=2 fi+1 

is a complex tensor of order R and size F x M 2 x • • • Mr. 
Similar expressions are evoked in, among others, [11]. The 
new tensor S can also be written as the concatenation of F 
tensors along the first dimension 

<S = <SiUi<S2LJi---Ui<Sf (28) 

where each Sf, f = 1,... ,F is a modal (i? — 1)-D signal of 
size 1 X M 2 X • • • X Mr containing a single (i? — 1)-D tone: 

Sf = Cf SLf,2®SLf^Z®- ■ -^SLf^R. (29) 

The singular value decomposition (SVD) of Y(i) yields 

Y(i) = USV" (30) 

where matrices U and V contain respectively the left and 

right singular vectors of Y(i), and S is a diagonal matrix 
containing the singular values ai,i = 1,..., min{Mi, M{} 
sorted in a decreasing order. As the number of components in 
[y is equal to F, then an approximation of Y(i), denoted by 
Y(i), can be obtained using the first F principal components 
of the SVD: 

Y(i)=UfSfVH (31) 

where U f (resp. V f) stands for the matrix formed with the 
first F columns of U (resp. V) and Df = diag((Ji ,... ,aF)- 
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It can be established from (8) and (31) that Ai and Uj’ 
span the same subspace, and thus there exists an unknown 
nonsingular matrix T that satisfies 

Ai = UfT. (32) 

Denote by M (resp. M) the matrix obtained from M by delet¬ 
ing the first (resp. last) row. By harnessing the Vandermonde 
structure of Ai, there exists a diagonal matrix D such that 
Aj = AiD. Since A^ = U^T and Ai = U^T, then 
UpT = Uj’TD, which proves that matrix T can be estimated 
by the eigenvectors of U^Uf- 

Thereby S can be estimated from the noisy data and Ai 
using equation (26) as follows 

S^y*A\, (33) 

then Sf,f = are extracted from S according 

to (28). Each yf = c/ay i ®---®Sif u can be estimated by 
y^f'^ = Sf»Bif^i. The sparse multigrid algorithm for single 

tone (STSM) can be applied on each y^P, / = 1,..., F to 
estimate the parameters of modes. However, we propose in 
the following to improve the separated components using an 
iterative technique. 


B. Improving the Estimation Accuracy 

It is clear from (33) that, in the noisy case, the error in 
estimating S (due to the estimation of Ai) will propagate 
when estimating the parameters a/_ 2 , • • ■, a/,K- Hence, we pro¬ 
pose to improve iteratively the mode estimates. The following 
procedure is executed to update estimates at each iteration 
7 = 0,..., A 

1) apply STSM to estimate a/ 2 , ■ • ■, ^f,R, / = 1, • ■ •) F 

{a/, 2 ,..., a/,fl} = STSM(yf, p,, r = 2,..., i?) 

(34) 


2) estimate c/a/,i,/ = l,...,Fby least squares using the 
already estimated a/,2, • ■ •, a/,77, / = 1, ..., F 

cTS^i = ((a/,fi K K a/,2)T)^ ( 35 ) 

3 ) compute y^ 

- (i) __- 

y^ = c/a/,1 (g)a/,2 0 • • • (8)a/,i7 ( 36 ) 


- fj') ^ {*~1) 

where yy = y} 


>77 (*) '77 




yP,f = i,...,F, nP and nP = y- 

f ■ This iterative scheme will be analyzed in the next 


Finally, the algorithm we propose (MTSM: Multiple Tones 
Sparse Method) is summarized in Algorithm 3. Note that no 
association step of R-D modes is required. The initialization 
step consists in initializing: i) Ai and S using (32), (33) 
and (28), ii) the estimated single tones yP, f = 1,..., F. 
Note that the columns of Ai are iteratively updated without 
extracting the related modes, whereas the modes of the other 
dimensions are extracted at each iteration using (34). Solely 
after the last iteration (i = K), the parameters of the hrst 


Algorithm 3: Sparse multigrid method for multiple tones 
estimation (MTSM) 

input : A tensor y S ,? 7 q,) g N X N 

output: Parameters of the multiple i?-D modes : 

initialization: 

1) Compute Ai and <S/, / = l,...,i^ using (32), (33) and (28) 

2) yf=.S/.a/,i,/=l,...,F 

For / = 1,. .., F, compute yP using (34), (35) and (36) 

for i = 1 : iC do 

for / = 1 : F do 

yf = yp^'’ +'R-Pi 

compute using (34), (35) and (36) 

TZ.y = yf - yf, if / = F, then TZ-f 

end 

end 

For / = 1,..., F, extract a/, j using 

a/,1 = STSM{yp^ +'R.p\ri^,ric„r = 1) 
return {d/,^}/^^^^! 


dimension are extracted using STSM algorithm. K denotes 
the maximum number of iterations, which is fixed to 2 in the 
simulations since no improvement was observed for K > 2. 

C. Analysis of the Algorithm 

Following the separation step described in (30)-(33), we 
can state that the algorithm yields the expected solution when 
the SNR is sufficiently high. We want now to prove that 
the second stage (next iterations), in addition to estimating 
the parameters from the single tones, is also improving the 
estimation accuracy. The general idea is inspired from greedy 
forward/backward sparse approximation, where the solution is 
refined by adding/removing atoms to/from the set of activated 
atoms. The improvement of the estimates is stated by the 
following theorem. 

Theorem 3: Assuming that the noise S is sufficiently small 
such that the ordering of the singular values in S in (30) is 
the same as the ordering of the corresponding singular values 
when S — 0. Using the procedure expressed by (34), (35) and 
(36) to estimate y f at iteration i = 0,..., K 

yP = argmin - X\\ (37) 

XGH 

where R = {X G ((^Mi x-'-xMi? _ bi 0 b 2 0 • • • 0 b/j, 
hr G V for r y 1} with V = {v G C'^''|v = 
[1, u,..., V = exp(/3 -I- juj), /3 e R”, w G [0, 27r)}. 

Then, at each iteration i,i = the residual is de¬ 

creased: 

~ ^ (i) ~ ^ (i— 1) 

y-y < y-y ( 38 ) 

^ (■i) _ fp - (i) 

where ^ ^/ ■ 

Proof See Appendix B. ■ 

Figure 3 depicts a comparison between results obtained by the 
MTSM algorithm with two different values of A G (0, 2}. The 
results show that MTSM with the improving step yields accu¬ 
rate estimates as compared to MTSM without the improving 
step. 



for i = 1,..., M, where 



-10 -5 0 5 10 15 20 

SNR (dB) 


Fig. 3. Frequency total root-mean square error for a 3-D signal containing 
3 modes with identical modes in two dimensions and close modes in the 
first dimension (Signal #4 in Table III). N(0) = 20,rj,^ = 21, = 

11, (Ml, M 2 , Ms) = (10,10,10). 100 Monte-Carlo. 

D. Identifiability 

Based on the assumptions under which Algorithm 3 is 
operating, the identifiability condition can be stated as F < 
Ml and min{M 2 ,M/j} > 2. In [33], the condition is 
andF< Oti \^]- 

We note that, when Mr > 4,r = the number 

of identifiable modes is slightly smaller than in [33], but the 
proposed algorithm is able to outperform the conventional 
methods in terms of computational complexity and accuracy. 
In addition, another advantage of the proposed algorithm is 
clear when the number of samples in one or more dimensions 
is less than 4 (i.e. Mr < 4), where identifiability in [33] is not 
satisfied. This latter case (i.e. 3r, Mr < 4) can be encountered 
in signal processing applications when the size of one or more 
diversities (dimensions in our formulation problem) is less 
than 4. 

VI. Cramer-Rao Lower Bounds for i?-D Cisoids in 

Noise 

In this section, we derive the expressions of the CRLB for 
the parameters of R-D damped exponentials in Gaussian white 
noise. We then give the CRLB in the cases of single damped 
and undamped R-D cisoids. We consider the R-D sinusoidal 
model given in (3). Let 

® = [wi,l • • • 102,1 ■■ ■U:f,R Otl,l ■ ■ ■Oil,R 

<22,1 • • ■ 01f,R Ai . . . Af 4>1 ■ ■ ■ 4>fY 

be the unknown parameter vector. The aim here is to derive 
the CRLB of the parameters in 0. 

The joint probability density function (pdf) of y is 

(39) 

where /x(0) is the noise-free part of y and 

y =[y(l, ■ • ■ , 1, 1), ■ • ■ , y(l, • ■ • , 1, Mr), 
yil,...,2,l),...,y{l,...,2,MF), 
...,y{Mi,...,MF)V. (40) 

The ith entry of fJ,{9) can be written as: 

F R 

/=i <'=1 


Tltr+lM, 


mod Mr, 


and [ J is the floor function. In the following, we derive the 
expressions of the CRLB in the general case (F > 1) and then 
we deduce the result corresponding to a single R-D modal 
signal (F = 1). 

A. Derivation of the CRLB 

Given the joint pdf in (39), the {k, 1) entry of the Fisher 
information matrix is [34], [35]: 




[F(0)]fci = ^Re 


We now express the derivatives diJ,{6)i/d0k for i = 1,..., M 
and k = 1, •.., 2RF + 2F. 

• For k = 1,..., RF, we have 






with r{k) = [{k — 1) mod i?] + 1 and f{k) = [{k — 
1)/R\ + 1. 

For k = RF + 1,..., 2RF: 


dn{e)i 


— ti,r(k)Cf(k) ' 


with r{k) = [{k — RF — 1) mod i?] + 1 and f{k) = 
[{k- RF -1)/R\ + 1. 

. For fc = 2RF + 1,..., 2RF + F: 

ti,r 

= ® ^46) 

r=l 

where ffk) = k — 2RF. 

. For fc = 2RF + F+1,..., 2RF + 2F: 


dii{9)^ 


~ n' 


where /(fc) = k — 2RF — F. 

Hence, the M x (2RF + 2F) matrix dii{9)/d9 expresses as 
= Z'^ Zcf, jZ</.].blkdiag(A,A,lF,A) 


Z' = [z;,..., Z'p] G with z'/(z, I) = u,i n a}-; 


A = blkdiag(AilF,..., AfIa) G 
$ = blkdiagje^'^^Ifi,..., G 

R 

Z = [zi,...,zf] G C^^-^.with z/(z) = 

r—1 

X = diag([Ai,..., Af]) G 
cj) = diag([e^‘^i,..., G 
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Finally, the inverse of the Fisher information matrix is 
F~\6») = yS-i [Re{V^V}]”^S-i = 

(55) 

where Re{-} stands for the real part. The CRLB of 9k is given 
by [F~^(0)]fcfc. More explicitly, for / = and r = 

2(7^Wi 

CRLB(a;/,0 = - ^5^^ 


CRLB(a/,^) = 


2cr2w 


'it(/-l)+r,it(/-l)+r 

X2 


CRLB(A/) = 2a‘^'W2RF+f,2RF+f 

2a'^'W2RF+F+f,2RF+F+f 


CRLB(()i/) = 




(57) 

(58) 

(59) 


Theorem 4: For the general i?-D exponential process, the 
CRLB’s for f = 1,..., F and r = 1,... ,R satisfy 

CRLB(a;/,^) = CRLB(a/,^) (60) 

CRLB(A/) = A2CRLB((/)/) (61) 

Proof It is based on the special block structure of matrix 
Re{V*^V} (see for instance [34]). 


B. Single Mode Case 

In this section, the CRLB’s will be simplified in the case 
of a single R-D modal signal {F = 1) to obtain more 
precise details on their parameter dependency. For the sake 
of simplicity, the subscripts denoting the mode / = 1 will 
be omitted. First, assume that |ar| = exp(ar) < 1. We 
shall express the products 7/^7J, and After some 

calculations, we get: 

[7,'^7.'u= n 


1 - la, 


\2Mr 


r—1 
r^n^k 

{ Mn-l 


1 - |a,.|2 

Mk-l 


X < 


E 

m—O 

Mn-1 

E 


1 2m 


E 

m—O 


m\ak\ 




V m—0 


z'^z ^ 


r—1 

R 


1 - la, 


|2Mr 


[z'''z]„=n 


r—1 

r^n 


1 - \ar\^ 
1 - 


1 — ttr- 


Mn-1 

X ^ TO|a„| 
m=0 


if k 

if n = k 

(62) 

(63) 

(64) 


Denoting = nf=i(l “ 


^M„-l 




la, 


|2Mr 
12m 


|arl ’')/(! |ar| )? 


^Mn-1 I 
jrn—0 I 

P™, we then obtain: 


and q 2 {n) = 


[P]„fc = M(“) X 


qi{n)qi{k), if n ^ k 
q 2 {n), if n = k 


G = 

[Q]„ = M(“)<?i(n), 


(65) 

( 66 ) 
(67) 


and 


Re{V”V} 


P 0 0 Q 

0 P Q 0 
0 QT G 0 
QT 0 0 G 


( 68 ) 


The inversion of RejV'^V} yields the following expressions 
of the CRLB’s: 


2 

CRLB(.9.) = CRLB(a.) = 

(i-ia,n^(i-Kr-)^ 

[-M2|a,|2M.(l _ |a,|2)2 + |a,|2(l _ |a,|2M.)2] ’ 

(69) 


CRLB(A) 


^ = cklb(^) = E(^ 


1 +E 


<lKr) 


^ 92(r) - qf{r) 


■ (70) 


Finally, for a single R-D purely harmonic signal (a^ = 0,Vr), 
we have = n^=i = M and taking the limit of the 

CRLB’s when ar ^ 0 leads to: 


is. CRI,BK)= y 


lim 

>-0 


CRLB (A) 


A2 


2X^M 


1 + 3E 


Mr - 1 


Mr 


1 


(71) 

(72) 


Hence, for the undamped case, our result in (71) is consistent 
with [11]. 


VII. Simulation Results 

Numerical simulations have been carried out to assess 
the performances of the proposed method for 2-D and 3- 
D modal signals in the presence of white Gaussian noise. 
The performances are measured by the total root-mean square 
error (RMSE) on estimated parameters and the computa¬ 
tional time. The total RMSE is defined as RMSEtotai = 

\lm'^P {Ef=i E/=i(Ot - 5/t)^} where is an esti¬ 
mate of ^f^r, and Ep is the average on p Monte-Carlo trials. In 
our simulations, ^f^r can be either a frequency or a damping 
factor. 

A. RMSE for 2-D and 3-D Signals 

Experiment 1: to show the interest of the multigrid scheme, 
this experiment presents the results obtained on Signal #1 with 
different multigrid levels and different initial grids. Signal #1 
is a single tone 2-D modal signal of size 10 x 10 whose 
parameters are presented in Table I. The number of multigrid 
levels is fixed to L = 2, i.e., £ = 0,1,2. Then the results are 
presented as a function of the number of atoms in the initial 
dictionaries N{0) and the number of atoms rj^, or pa added 
at each level £. The results we obtain for the first step, i.e., 
for the harmonic estimation, are presented in Eigure 4. We 
can observe that the frequency RMSE obtained with the R-D 
sparse algorithm can reach the CRLB using a uniform initial 
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TABLE I 

2-D PARAMETERS OF SIGNAL #1 


1 f 1 

1 J^f.l 1 

“/.l 1 

1 ’^f,2 1 

“.f,2 1 

1 Gf 1 

1 ^ 1 

1 0.22 1 

-0.011 1 

1 0.34 1 

-0.015 1 



harmonic dictionary of 10 atoms and rj,^ = 31 (Figure 4.a). 
Figure 4.b shows that the frequency RMSE is improved at low 
SNR if the initial dictionary contains 20 atoms, and reaches the 
optimal estimates with = 21. Figure 5 shows the damping 
factor RMSE obtained by R-D sparse algorithm using different 
settings of the initial damping factor dictionary and rj^. We 
can observe that the damping factor RMSE depends on the 
number of atoms in the dictionary, the more atoms the better. 
At low SNR, the RMSE also depends /^min- Therefore, it is 
better to choose /3min with small absolute value if we have 
a prior knowledge of the interval of damping factors in the 


signal. In general, the estimation error is of order Eor 

instance, in the frequency step estimation, we recommend to 

chose N{0) to be greater than or equal to if we want 

a good accuracy at lower SNR levels. Otherwise, we can set 

A^(0) = Mr- Once N{Q) is set, rj can be chosen with respect 

to the desired accuracy. Let e be the desired estimation error, 

then e = , and we can set n = , ^ 

I y/Tmfi) 

In the rest of this section, the proposed algorithms are com¬ 
pared with 2-D ESPRIT [ 6 ], Tensor-ESPRIT [10], PUMA [12] 
and TPUMA [11]. If the R-D signal contains one tone then 
Algorithm 2 (STSM) is used, otherwise Algorithm 3 (MTSM) 
is used. Thus, to facilitate notation, both proposed algorithms. 
Algorithm 2 and Algorithm 3, will be called R-D sparse. 
For the proposed method, the initial grid used to build the 
harmonic dictionary is the same for all dimensions; it contains 
50 frequency points uniformly distributed over the interval 
[0,1) and 10 damping factors (3 G [—0.05,0]. To simulate a 
random dictionary, at each run, the frequency grid is perturbed 
by a small random quantity. As a consequence of experiment 1, 
we use the following settings (L,p^,? 7 ^) = (2,21,11). The 
number of iterations in Algorithm 3 is set to AT = 2 because 
no improvement was observed for K > 2. 

Since the proposed method is applied directly on data 
without using spatial smoothing, i.e., it does not require the 
construction of a large matrix or an augmented order tensor, 
then a relevant comparison will be with algorithms that do 
not use spatial smoothing. Thereby, in the next experiments, 
the proposed algorithm is compared to PUMA [12] and 
TPUMA [11], which are algorithms that do not require spatial 
smoothing. We also report comparisons with 2-D ESPRIT [ 6 ] 
and Tensor-ESPRIT [10], which need spatial smoothing. 


1) Single tone R-D modal signal: 

Experiment 2: This experiment tends to show the efficiency 
of the proposed algorithm in estimating parameters of single 
tone R-D modal signals. We simulate a 2-D signal of size 
10 X 10 (Signal #1) whose parameters are presented in Table I. 
Our R-D sparse algorithm is compared to 2-D ESPRIT [ 6 ] and 
PUMA [12]. Eor each level of noise, 1000 Monte-Carlo trials 
are performed. Eigure 6 shows the obtained results. We can 
observe that: i) the proposed algorithm and PUMA reach the 


10 


I 

LU 10'' 
03 


10 ' 
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O N{0) = 10, rii, = 11 
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- - CRLB 
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(a) A^(0) = 10 



* iV(0) = 20, Vu = ^ 

0 N(0) = 20, riu = 11 
n Ar(0) = 20, riu = 21 

-CRLB 




(f) 

tr 
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SNR{dB) 


(b) 7V(0) = 20 

Fig. 4. Frequency RMSE using i?-D sparse algorithm with different rji'- 2- 
D signal containing a single tone (Signal #1). (Mi, M 2 ) = (10, 10). 1000 
Monte-Carlo trials, (a) The initial harmonic dictionary contains N(0) = 10 
atoms, (b) N{0) = 20 atoms. 


TABLE II 

Different configurations for experiments 3 through 5 



F 

Dimension 1 

Dimension 2 

Dimension 3 

Exp. 3 

2 

> Apr 

Ah' > Apr 

Ah' > Apr 

Exp. 4 

3 

Au > Apr 

3 identical modes 

3 identical modes 

Exp. 5 

3 

Au <C Apr 

3 identical modes 

3 identical modes 


CRLB and outperform 2-D ESPRIT, ii) R-D sparse outperform 
PUMA in SNR less than 3 dB. 

2) Multiple tones R-D modal signals: Several configu¬ 
rations are studied in the case of multiple tones to com¬ 
pare the proposed algorithm with Tensor-ESPRIT [10] and 
TPUMA [11]. These configurations (Experiments 3, 4, 5) 
are summarized in Table II, in which the number of modes 
and the distance between frequencies in different dimensions 
are varied. Apr denotes the Rayleigh frequency resolution 
limit, which has the same value in all dimensions because 
Ml = M 2 = M 3 . In Experiment 6 we examine the case 
when the size of only one dimension is larger than 4, i.e., the 
identifiability condition of [33] is not satisfied. The parameters 
of the used signals are given in Table III. 

Experiment 3: In this experiment, we simulate a 3-D signal 
(Signal #2) of size 8 x 8 x 8 and containing two modes whose 
frequencies in each dimension are well separated. Parameters 
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TABLE III 

3-D PARAMETERS OF SIGNAL #2 THROUGH #5 


Signal 

^.f.l 

“f,i 

''f.2 

“.f.2 


a/,3 


#2 

0.40 

-0.01 

0.1 

-0.01 

0.1 

-0.01 

1 


0.20 

-0.01 

0.3 

-0.15 

0.25 

-0.01 

1 

#3 

0.30 

-0.01 

0.31 

-0.01 

0.22 

-0.01 

I 


0.10 

-0.01 

0.45 

-0.015 

0.11 

-0.01 

1 


0.20 

-0.01 

0.31 

-0.01 

0.11 

-0.01 

I 

#4 

0.28 

-0.01 

0.31 

-0.01 

0.22 

-0.01 

I 


0.12 

-0.01 

0.45 

-0.015 

0.11 

-0.01 

1 


0.20 

-0.01 

0.31 

-0.01 

0.11 

-0.01 

1 

#5 

0.30 

-0.01 

0.1 

-0.01 

0.1 

-0.01 

I 


0.13 

-0.01 

0.45 

-0.015 

0.4 

-0.01 

1 


0.20 

-0.01 

0.31 

-0.01 

0.1 

-0.01 

1 


0.42 

-0.012 

0.22 

-0.01 

0.32 

-0.01 

1 


Fig. 5. Damping factor RMSE using R-D sparse algorithm with different r/c 
/3niin and Ai(0). 2-D signal containing a single tone (Signal#!). {Mi,M 2 ) = 
(10,10). 1000 Monte-Carlo trials. 
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Fig. 6. Frequency total root-mean square error for a 2-D signal containing a 
single tone (Signal #1). {Mi,M 2 ) = (10,10). 1000 Monte-Carlo trials. 

of the signal are given in Table III. Figure 7 shows the obtained 
results. Here, the proposed method performs as TPUMA. 
Tensor-ESPRIT yields slightly worse RMSE. 

Experiment 4: 3-D signal of size 10 x 10 x 10 containing 
three 3-D modes (Signal #3). Note that there exists identical 
modes in two dimensions and frequencies in the first dimen¬ 
sion are separated by 1/Mi. Eigure 8 shows the results. In 
this configuration, TPUMA and Tensor-ESPRIT give similar 
results and the proposed method performs better for all SNR 
levels. 

Experiment 5: 3-D signal of size 10 x 10 x 10 containing 
three 3-D modes (Signal #4). Note that there exists iden¬ 
tical modes in two dimensions and frequencies in the first 
dimension are separated by less than 1/Mi. The results are 
shown on Eigure 9. Here again the proposed R-D sparse 
approach performs better than TPUMA and Tensor-ESPRIT. 
Observe also that Tensor-ESPRIT outperforms TPUMA in 
this configuration (close frequencies and identical modes in 
dimensions 2-3). 

Experiment 6: Results on Signal #5 of size 10 x 3 x 3 
containing 4 modes are given in Eigure 10. We observe that 
the proposed method outperforms TPUMA algorithm mainly 
in low SNR levels. 


B. Numerical Complexity 

It is known that in the case of 1-D signals of size M, 
OMP costs 0{NFM) in terms of multiplications [36]; F 
is the sparsity (number of components) and N is the num¬ 
ber of atoms in the dictionary. Eor a M-measurements R- 
D signal, the complexity of the STSM algorithm over a 
set of L multigrid levels is 0{MNLR), assuming that the 
number of dictionary atoms is maintained constant (equal 
to N) over all levels. Regarding the approach proposed in 
Algorithm 3, the main operations are the call of STSM and the 
update of c/a/d = '^/(i) Kl • • • Kl which has a 

complexity of 0 (M) since {{&f,R Kl • • • Kl ay 2 )^)^ is a row of 
length n ^2 and is a matrix of size Mi x 0^=2 ^^r- 
Therefore, the whole complexity of the proposed algorithm is 
O {{NL{F(^R — 1)K -f 1) -f FK)M), which is linear in the 
number of measurements M. The complexity of the Tensor- 
ESPRIT algorithm with spatial smoothing is mainly related to 
that of the SVD which is at least 0{ktF{R + 1)PM) where 
kt is a constant depending on the implementation of the SVD 
algorithm. Here P = 0^=1 Pr where {Pr-}^=i are design 
parameters used to get smoothed measurements (see [ 10 ]). 
The accuracy of the estimates provided by ESPRIT depends on 
these parameters. Since the optimal value for Pr is a fraction of 
Mr (e.g. [37]-[39]), the complexity of the SVD step is, in fact, 
O(M^). The complexities of PUMA and TPUMA algorithms 
are 0{M^) and 0{ktM{R+F-l))+J2^^^ 0{K{F+1)M^), 
respectively. Compared to PUMA and TPUMA, the proposed 
algorithm has an attractive computational complexity for large 
size signals. Eigure 11 shows the CPU time results of the 
proposed, Tensor-ESPRIT and TPUMA algorithms versus 
Ml for a 3-D damped signal containing two modes with 
M 2 = M 3 = 4. We observe that the proposed method has low 
computational complexity compared to TPUMA and Tensor- 
ESPRIT when Ml is large. 

VHI. Conclusion 

We presented an efficient sparse estimation approach for 
the analysis of multidimensional (R-D) damped or undamped 
modal signals. The idea consists in exploiting the simultaneous 
sparse approximation principle to separate this joint estimation 
problem into R multiple measurements problems. To be able 
to handle large size signals and yield accurate estimates, a 
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Fig. 7. Frequency total root-mean square error for a 3-D signal containing Fig. 9. Frequency total root-mean square eiTor for a 3-D signal containing 3 
two 3-D modes (Signal #2). (Mi, M 2 , M^) = (8,8,8). 1000 Monte-Carlo. modes with identical modes in two dimensions (Signal #4), same as Signal 

#4 with close modes in the first dimension (0.28,0.12,0.2). (Mi, M 2 , M3) = 
(10,10,10). 200 Monte-Carlo. 
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multigrid dictionary refinement scheme is associated with the 
simultaneous orthogonal matching pursuit (SOMP) algorithm. 
We gave the convergence proof of the the rehnement procedure 
in the single tone case. Then, for the general multiple tones 
R-D case, the signal tensor model is decomposed in order 
to handle each tone separately in an iterative scheme so that 
the pairing of the R-D parameters is automatically achieved. 
Also, the CRLB of the i?-D modal signal parameters were 
derived. The tests performed on simulated signals showed that 
the proposed algorithm attains the CRLB and outperforms 
state-of-the-art subspace algorithms. We also have shown that 
the complexity of the algorithm is linear with respect to the 
number of measurements, which allows the processing of large 
size signals. Finally, it is worth mentioning that this approach 
can be straightforwardly applied to other multidimensional 
array processing problems. 

Appendix A 
SOMP Algorithm 

In this appendix we report the SOMP method (Algo¬ 
rithm 4) [23]. In this description, jmj denotes the TO 2 * vector 
of the canonical basis in . 


Algorithm 4: SOMP 

input : A matrix Y G matrix Q G (with 

normalized columns) 

output: An index set Q of activated atoms. A matrix of sparse vectors 

X G C^x-^2 

initialization: = 0, = 0, X = 0, Rq = Y 

while halting criterion false do 
k = k 1 

nk G arg max„J]m2=l I (R-fc-lJma . Hr.) | 
fl/j — U {rife} 

Xfc = (Q« Y 

R-fc = Y — Qn^Xfc 

end 

return O = Ofc, X = Xj, 


Appendix B 
Proof of Theorem 3 

We begin the proof by introducing the following lemma. 
Lemma 1: Consider y = V-f-AV, where y is the perturbed 
version of the data tensor y and Ay is the perturbation. 
Assuming that AV is sufficiently small such that the ordering 
of the F singular values in S in (30) is the same as the 
ordering of the corresponding singular values when Ay = 0. 
Then the perturbation Ay f contains a linear combination of 
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Fig. 11. Average CPU time for a single run under M 2 = M3 = 4 and 
F = 2. 


all yfJ = l,...,F: 


F 

Ay/f = T>f + '^vf^^y^ 
2=1 


where vj = [vf,i,...,Vf,F] = AA|(/, ;)Ai and X>/ = 
Ay,eif,,A\{f,:)+ysj*Aaf,,. 

Proof From (33) S = 3^#A|, we differentiate and obtain 

AS = Ay*A\ aa| 

1 1 

Then, 


ASf = Ay. Alif,:)+S. AAl (/,:) Ai 
1 1 ^ ^ ^ 

VT 

= Ay.Al{f,:)+S.y 

F 

^ Ay .Alif,:) +Y^ Vf,pSp 

p=i 


Therefore, a/ 2 ) • ■ • j &/,«> / = can be estimated 

using STSM algorithm since 

= c/a/_r(a/,fl Kl • • • Kl a.f^r+1 a/,r-il^ 

■■■^(a/,i+P/,/a/.i + Aa/,i)) 



A>;.a;.iAl(/,:) 


(r) 


Since has the following form 


+ ^/./a/,1 + Aa/,i)(a/,^j K K a/, 2 ) 


^ u/.p^p + Ay.a/.iAU/,:) 

I p=i,p// 


( 1 ) 


we estimate c/a/,i by least squares once a/, 2 , • ■ •, a/,/j are 
estimated using STSM 


c/a^ = min ||y/°^ — a® a /,2 


5%, it I 


So, we put yj = c/a /,1 (g) a /,2 ( 8 ) • • • ( 8 ) and TZf = 
- (0) ~ (0) 

yr -Ff . Therefore, the procedure to estimate i^/ at 
iteration i = 0,..., K can be summarized in (34), (35) and 
(36). Note that this procedure is optimal because STSM and 
the least squares are optimal when they are used to estimate 
a/, 2 ,..., a/,/{,/ = 1 ,...,F and c/a/,i,/ = 1 ,...,F, 

respectively. 

Now we present the technique for improving the estimation 
of yf. Let nf = = y - E/=1 yf and 

y^P = argmin \\yf^ + - A|| (73) 

A-ew 

where = yf + TZfl^ - yf, f = 1,..., F, and yf 

is an improved estimate of y f. We follow the same procedure 

~ ( 1 ) 

as described in equations (34), (35) and (36) to calculate ^/ . 
We can state that there is improvement in the estimation of 

yf if 


yf is estimated using yf = <S/*a/,i, we differentiate and 
obtain 


y-Ey 


( 1 ) 


/=! 


< 


y-T,yi 

/=! 


(0) 


(74) 


A3^/ = ASf • a /,1 + 5/ • Aa /,1 

F 

= H ^f’pyp + *^7 T +^y* a/.iAl(/, 0 


Using the previous lemma 



*5/ '(a/,! + ti/,/a/,i + Aa/) + 

F 

^ u/,p3^p + A>^.a/,iAl(/,:) 

p=i.p#/ 


We have ||7ei^)|| = - yf + E/= 2 ^/ + V|| where 

V = V — y and y = Fifiyf ■ ‘-an be verified that 


||7?.W|| = 


yf+EW-yr)+ E ^p+v|-yy 

i p=i p=/+i / 


( 1 ) 


117 ?.' 




yf + E(yf - yf) + E ^p+V1 - y 


,(i) 


p=i 


p=/+i 


,(i) ■ 


However, from equation (73), 3^/ is the minimizer with 
respect to AT £ 7f of 

Eiyr - S'?’) + E f - ■» 

P=1 p=/+l / 
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Therefore, < ||7?.^^2ill) / = 1,...,-F. As con¬ 
sequence, ||7?.y‘^|| < ||7?.^^||, which we are seeking in 

expression (74). Similarly, we can prove that ||7?.^^|| < 

11II, i > 1, using the general forms of and 


Tt 


- ( y? + - yf) + E +v 

\ p=i p=/+i ) 

-yf 


/-I 


p=i 


,( 0 ) _ 


which we are seeking in (38). 



(75) 

+ E (S’?’ 

p=/+i 

-vr'Vv 


(76) 
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